#pragma once

// 数值计算的精度最好取 double，不然精度不够
using Real = double;

#include <iostream>
#include <initializer_list>

namespace LinearAlgebra
{
    // 预先声明
    class Vector;
    class Matrix;

    /**
     * @file Matrix.h
     * @author xingyifan
     * @date 2022-11-28 21:51
     *
     * @description: 向量类型
     */

    class Vector
    {
    protected:
        const int m_dim;
        Real *m_value;

    public:
        /**
         * @file Matrix.h
         * @author xingyifan
         * @date 2022-11-28 21:54
         *
         * @description: 构造函数
         */

        Vector(int D, Real t = 0);             // 初始化构造函数
        Vector(std::initializer_list<Real> l); // 可以通过列表初始化
        Vector(const Vector &rhs);             // 拷贝构造

        /**
         * @file Matrix.h
         * @author xingyifan
         * @date 2022-11-30 09:25
         *
         * @description: 向量信息
         */

        int size() const; // 获得向量长度

        /**
         * @file Matrix.h
         * @author xingyifan
         * @date 2022-11-28 21:57
         *
         * @description: 操作符重载
         */

        Real &operator[](int i);      // 返回引用
        Real operator[](int i) const; // 返回值

        // 向量加减
        Vector operator+(const Vector &rhs) const;
        Vector operator-(const Vector &rhs) const;

        // 常数乘除法
        Vector operator*(Real t) const;
        Vector operator/(Real t) const;

        // 向量乘矩阵
        Matrix operator*(const Matrix &rhs) const;

        // 赋值
        void operator=(const Vector &rhs);

        // 取反操作
        Vector operator-() const;

        // 输出函数，此函数就放在这里
        friend std::ostream &operator<<(std::ostream &o, const Vector &v)
        {
            for (int i = 0; i < v.size(); i++)
                o << v[i] << std::endl;
            return o;
        }

        /**
         * @file Matrix.h
         * @author xingyifan
         * @date 2022-11-28 22:12
         *
         * @description: 基本向量运算
         */

        void sort();                       // 冒泡排序
        Vector reverse();                  // 反转向量
        int find(Real t, int i = 0) const; // 从 i 开始查找元素，返回第一个下标
        Real sum() const;                  // 元素求和
        Vector abs() const;                // 取绝对值
        Vector sgn() const;                // 获取符号向量
        Real norm(int nt = 2) const;       // 向量范数
        Real minimum() const;              // 取最小元
        Real maximum() const;              // 取最大元
        int maxIndex() const;              // 取具有最大模的分量
        Matrix toMatrix() const;           // 转化成矩阵
        Matrix transpose() const;          // 转置成矩阵

        /**
         * @file Matrix.h
         * @author xingyifan
         * @date 2022-11-30 09:27
         *
         * @description: 释放资源
         */

        ~Vector();
    };

    /**
     * @file Matrix.h
     * @author xingyifan
     * @date 2022-11-28 22:29
     *
     * @description: 矩阵类型
     */

    class Matrix
    {
    protected:
        const int m_row, m_col;
        Real **m_value;

    public:
        /**
         * @file Matrix.h
         * @author xingyifan
         * @date 2022-11-28 22:30
         *
         * @description: 构造函数
         */

        Matrix(int Row, int Col, Real t = 0);                    // 初始化构造函数
        Matrix(int Row, int Col, std::initializer_list<Real> l); // 通过列表初始化
        Matrix(const Matrix &rhs);                               // 拷贝构造

        /**
         * @file Matrix.h
         * @author xingyifan
         * @date 2022-11-30 09:52
         *
         * @description: 矩阵信息
         */

        int row() const; // 获得行数
        int col() const; // 获得列数

        /**
         * @file Matrix.h
         * @author xingyifan
         * @date 2022-11-28 22:37
         *
         * @description: 操作符重载
         */

        Vector operator()(int j) const;      // 返回第 j 列
        Real operator()(int i, int j) const; // 返回 i,j 元素
        Real &operator()(int i, int j);      // 返回 i,j 元素引用

        // 矩阵加减法，矩阵乘向量
        Matrix operator+(const Matrix &rhs) const;
        Matrix operator-(const Matrix &rhs) const;
        Matrix operator*(Real t) const;
        Matrix operator/(Real t) const;
        Vector operator*(const Vector &v) const;
        Matrix operator*(const Matrix &rhs) const;

        void operator=(const Matrix &rhs); // 直接赋值
        Matrix operator-() const;          // 取反操作

        // 输出函数，此函数就放在这里
        friend std::ostream &operator<<(std::ostream &o, const Matrix &m)
        {
            for (int i = 0; i < m.row(); i++)
            {
                for (int j = 0; j < m.col(); j++)
                    o << m(i, j) << " ";
                o << std::endl;
            }
            return o;
        }

        /**
         * @file Matrix.h
         * @author xingyifan
         * @date 2022-11-28 23:05
         *
         * @description: 元素操作
         */

        void insert(int j, Vector &v); // 在 j 列插入列向量 v
        Real norm(int nt = 0) const;   // 矩阵范数
        Real norm(char c) const;       // e 范数 / f 范数
        Matrix upper() const;          // 获得上三角
        Matrix lower() const;          // 获得下三角
        Matrix transpose() const;      // 转置

        /**
         * @file Matrix.h
         * @author xingyifan
         * @date 2022-11-30 10:34
         *
         * @description: 矩阵算法，不改变矩阵
         */

        Matrix inverse() const;                                            // 矩阵求逆
        Matrix LU() const;                                                 // 标准 LU 分解
        Matrix LU(int perm[]) const;                                       // 列主元 LU 分解，传入下标数组，返回置换后的数组
        Matrix Cholesky() const;                                           // Cholesky LL'
        Matrix CholeskyD() const;                                          // 优化的 Cholesky LDL'
        Matrix QR(Vector &d) const;                                        // QR 分解
        Matrix Hessenberg(Matrix *Q = nullptr) const;                      // 上 Hessenberg 化，在 Q 中累计正交变换
        Matrix Tridiagonal(Matrix *Q = nullptr) const;                     // 对称阵化为三对角对称阵，在 Q 中累计正交变换
        Matrix Didiagonal(Matrix *U = nullptr, Matrix *V = nullptr) const; // 矩阵上二对角化，在 U,V 中累计正交变换

        /**
         * @file Matrix.h
         * @author xingyifan
         * @date 2022-12-01 15:32
         *
         * @description: 矩阵迭代法，将对矩阵本身进行迭代
         */

        void HA(int startOfRow, int endOfRow, Vector &v);                                  // 计算指定范围的 Householder 变换，作用结果为 HA
        void AH(int startOfCol, int endOfCol, Vector &v);                                  // 计算指定范围的 Householder 变换，作用结果为 AH
        void HAHs(int startOfRow, int endOfRow, Vector &v);                                // 计算对称阵指定范围的 Householder 变换，作用结果为 HAH
        void GA(int i, int j, Vector &cs);                                                 // 计算指定分量的 Givens 变换，得到 GA
        void AG(int i, int j, Vector &cs);                                                 // 计算指定分量的 Givens 变换，得到 AG
        void AGT(int i, int j, Vector &cs);                                                // 计算指定分量的 Givens 变换，得到 AGT
        void doubleQR(int l = 0, int m = 0, Matrix *Q = nullptr);                          // 双重步 QR 迭代，l,m 是 H22 前面和后面分块的阶数
        void WilkinsonQR(int l = 0, int m = 0, Matrix *Q = nullptr);                       // Wilkinson QR 迭代，针对三对角对称矩阵
        void WilkinsonSVD(int l = 0, int m = 0, Matrix *U = nullptr, Matrix *V = nullptr); // Wilkinson 位移的 SVD 迭代

        /**
         * @file Matrix.h
         * @author xingyifan
         * @date 2022-11-30 09:40
         *
         * @description: 释放资源
         */

        ~Matrix();
    };

    /**
     * @file Matrix.h
     * @author xingyifan
     * @date 2022-11-29 09:39
     *
     * @description: 动态数组，用于稀疏矩阵
     */

    template <class T>
    class DynamicArray
    {
    public:
        /**
         * @file Matrix.h
         * @author xingyifan
         * @date 2022-11-29 09:39
         *
         * @description: 内置节点类型
         */

        struct Node
        {
            T value;
            Node *prev;
            Node *next;

            Node(const T &value = T(), Node *prev = nullptr, Node *next = nullptr) : prev(prev), next(next), value(value) {}
        };

    protected:
        /**
         * @file Matrix.h
         * @author xingyifan
         * @date 2022-11-29 09:40
         *
         * @description: 保有节点头部和数组长度
         */

        int m_size;
        Node *m_head;
        Node *m_tail;

    public:
        /**
         * @file Matrix.h
         * @author xingyifan
         * @date 2022-11-29 09:41
         *
         * @description: 构造函数
         */

        DynamicArray();                           // 直接构造
        DynamicArray(std::initializer_list<T> l); // 初始化构造
        DynamicArray(const DynamicArray &arr);    // 拷贝构造

        /**
         * @file Matrix.h
         * @author xingyifan
         * @date 2022-11-29 10:26
         *
         * @description: 获得数组信息
         */

        int size() const { return m_size; }
        bool empty() const { return (m_size == 0) ? true : false; }

        /**
         * @file Matrix.h
         * @author xingyifan
         * @date 2022-11-29 10:04
         *
         * @description: 操作符重载
         */

        T &operator[](int i);
        T operator[](int i) const;

        // 拷贝赋值
        void operator=(const DynamicArray &arr);

        /**
         * @file Matrix.h
         * @author xingyifan
         * @date 2022-11-29 10:09
         *
         * @description: 元素操作
         */

        T front() const;                   // 获得头部元素
        T back() const;                    // 获得尾部元素
        Node *node(int i) const;           // 获得指定的节点
        void insert(int i, T value);       // 在 i 位置插入元素
        Node *insert(Node *node, T value); // 在指定节点处插入元素并返回新节点
        void push_back(T value);           // 最后加入元素
        void push_front(T value);          // 最前面加入元素
        T pop_back();                      // 弹出最后一个元素
        T pop_front();                     // 弹出第一个元素
        int find(T value) const;           // 找到元素
        bool kill(T value);                // 删除指定元素
        T remove(int i);                   // 移除指定索引处的元素并返回
        void remove(Node *node);           // 移除指定节点的元素
        void clear();                      // 删除数组元素

        /**
         * @file Matrix.h
         * @author xingyifan
         * @date 2022-11-29 10:27
         *
         * @description: 释放内存
         */

        ~DynamicArray();
    };

    /**
     * @file Matrix.h
     * @author xingyifan
     * @date 2022-11-29 14:09
     *
     * @description: 稀疏矩阵类型
     */

    class SpMatrix
    {
    protected:
        const int m_row, m_col;       // 尺寸
        DynamicArray<int> m_IA, m_JA; // 存放索引
        DynamicArray<Real> m_AA;      // 对角元放在每行首位

    public:
        friend class LinearSolver; // 定义友元，方便访问

        /**
         * @file Matrix.h
         * @author xingyifan
         * @date 2022-11-29 14:10
         *
         * @description: 构造函数
         */

        SpMatrix(int Row, int Col);         // 直接构造
        SpMatrix(const SpMatrix &spm);      // 拷贝构造
        explicit SpMatrix(const Matrix &m); // 用普通矩阵显式构造，防止隐式转换

        /**
         * @file Matrix.h
         * @author xingyifan
         * @date 2022-11-29 15:45
         *
         * @description: 结构信息
         */

        int row() const;                   // 获得行数
        int col() const;                   // 获得列数
        int explicitNum() const;           // 获得显式表达元素的个数
        Real diag(int i) const;            // 获取对角元
        int startOfRow(int i) const;       // 获取指定行的开头
        int indexOfCol(int j) const;       // 通过全局指标获取元素所在的列位置
        Real global(int j) const;          // 通过全局指标获取元素值
        int globalCol(int i, int j) const; // 获取(i,j)元素的全局列标

        /**
         * @file Matrix.h
         * @author xingyifan
         * @date 2022-11-29 15:51
         *
         * @description: 插入和修改元素
         */

        Real operator()(int i, int j) const; // 获得存在的元素，不存在的元素认为是零
        Real &operator()(int i, int j);      // 获得元素引用，不存在的元素就创建元素

        SpMatrix compress() const;  // 压缩矩阵的零元素，返回压缩矩阵
        SpMatrix transpose() const; // 转置矩阵
        Matrix toMatrix() const;    // 转化为一般矩阵

        /**
         * @file Matrix.h
         * @author xingyifan
         * @date 2022-11-29 17:14
         *
         * @description: 操作符重载
         */

        // 加减法
        SpMatrix operator+(const SpMatrix &spm) const;
        SpMatrix operator-(const SpMatrix &spm) const;

        // 乘法运算
        SpMatrix operator*(Real k) const;
        SpMatrix operator/(Real k) const;

        // 稀疏矩阵乘稀疏矩阵
        SpMatrix operator*(const SpMatrix &spm) const;

        // 拷贝赋值
        void operator=(const SpMatrix &spm);

        // 输出函数，此函数就放在这里
        friend std::ostream &operator<<(std::ostream &o, const SpMatrix &A)
        {
            int n = A.row();

            // 使用节点加速访问
            DynamicArray<int>::Node *node_IA = A.m_IA.node(0);
            for (int i = 0; i < n; i++)
            {
                DynamicArray<int>::Node *node_JA = A.m_JA.node(node_IA->value);
                DynamicArray<Real>::Node *node_AA = A.m_AA.node(node_IA->value);

                // IA[m] = nnz，边界值
                for (int j = node_IA->value; j < node_IA->next->value; j++)
                {
                    o << i << ", " << node_JA->value << ", " << node_AA->value << std::endl;

                    node_JA = node_JA->next;
                    node_AA = node_AA->next;
                }
                node_IA = node_IA->next;
            }
            return o;
        }
    };

    /**
     * @file Matrix.h
     * @author xingyifan
     * @date 2022-11-30 10:28
     *
     * @description: 基本线性求解器
     */

    class LinearSolver
    {
    public:
        /**
         * @file Matrix.h
         * @author xingyifan
         * @date 2022-11-30 10:29
         *
         * @description: 标准求解线性系统
         */

        Vector solveUpper(const Matrix &A, const Vector &b, bool unit = false); // 是否按照单位上三角阵求解
        Vector solveLower(const Matrix &A, const Vector &b, bool unit = false); // 是否按照单位下三角阵求解
        Vector solve(const Matrix &A, const Vector &b, bool choose = false);    // 是否选主元求解

        /**
         * @file Matrix.h
         * @author xingyifan
         * @date 2022-11-30 13:34
         *
         * @description: 对称线性系统求解
         */

        Vector Cholesky(const Matrix &A, const Vector &b, bool improve = false); // 是否优化为 LDL 分解

        /**
         * @file Matrix.h
         * @author xingyifan
         * @date 2022-12-06 16:03
         *
         * @description: 矩阵信息
         */

        Real condition(const Matrix &A); // 计算条件数，基于无穷范数

        /**
         * @file Matrix.h
         * @author xingyifan
         * @date 2022-11-30 13:34
         *
         * @description: 古典迭代法
         */

        Vector Jacobi(const Matrix &A, const Vector &b, Real omega = 1); // Jacobi 迭代法
        Vector GaussSediel(const Matrix &A, const Vector &b);            // Gauss-Sediel 迭代
        Vector SOR(const Matrix &A, const Vector &b, Real w);            // 超松弛迭代

        /**
         * @file Matrix.h
         * @author xingyifan
         * @date 2022-11-30 13:44
         *
         * @description: 稀疏矩阵的迭代法
         */

        Vector Jacobi(const SpMatrix &A, const Vector &b, Real omega = 1); // Jacobi 迭代法
        Vector GaussSediel(const SpMatrix &A, const Vector &b);            // Gauss-Sediel 迭代
        Vector SOR(const SpMatrix &A, const Vector &b, Real w);            // 超松弛迭代

        /**
         * @file Matrix.h
         * @author xingyifan
         * @date 2022-11-30 10:31
         *
         * @description: 最小二乘问题
         */

        Vector QR(const Matrix &A, Vector &b);       // QR 分解法
        Vector CG(const Matrix &A, const Vector &b); // 共轭梯度法

        /**
         * @file Matrix.h
         * @author xingyifan
         * @date 2022-11-30 16:38
         *
         * @description: 特征值和特征向量
         */

        Real power(const Matrix &A, Real mu = 0);                                    // 幂法计算 A-muI 的最大特征值，然后加回 mu
        Vector inversePower(const Matrix &A, Real lambda);                           // 反幂法计算指定特征值的特征向量（半次迭代法）
        Matrix implicitQR(const Matrix &A, Matrix *Q = nullptr, bool ifSym = false); // 隐式 QR 方法

        /**
         * @file Matrix.h
         * @author xingyifan
         * @date 2022-12-06 16:10
         *
         * @description: 对称矩阵的特征值和特征向量
         */

        Matrix passJacobi(Matrix &A, Real sigma); // 过关 Jacobi 方法，将 A 变换为特征值对角阵，返回正交阵 Q，每一列是特征向量
        Vector bisection(const Matrix &A);        // 二分法，计算三对角对称阵的全部特征值

        /**
         * @file Matrix.h
         * @author xingyifan
         * @date 2022-12-06 20:18
         *
         * @description: 奇异值分解
         */

        Matrix SVD(const Matrix &A, Matrix *U = nullptr, Matrix *V = nullptr); // 返回分解后的对角阵，正交阵存放在 U,V 中
    };

    /**
     * @file Matrix.h
     * @author xingyifan
     * @date 2022-11-30 21:36
     *
     * @description: 基本向量操作和矩阵操作
     */

    Real dot(const Vector &vec1, const Vector &vec2); // 向量点积
    Real Householder(Vector &x, Vector &v);           // Householder 变换
    Vector Givens(Vector &v, int i = 0, int j = 1);   // Givens 变换
    Matrix Hilbert(int n);                            // 希尔伯特矩阵
    void rotate(Matrix &A, Matrix &Q, int p, int q);  // 对称矩阵绕 p,q 分量旋转 JT*A*J，Q 累计 J
};